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Abstract 

In  this  paper  we  discuss  the  numerical  integration  of  Lie-Poisson  Systems  using  the  mid¬ 
point  rule.  Since  such  systems  result  from  the  reduction  of  hamiltonian  systems  with 
symmetry  by  Lie  group  actions,  we  also  present  examples  of  reconstruction  rules  for  the 
full  dynamics.  A  primary  motivation  is  to  preserve  in  the  integration  process,  various 
conserved  quantities  of  the  original  dynamics.  A  main  result  of  this  paper  is  an  0  (/?3 )  error 
estimate  for  the  Lie-Poisson  structure  where  h  is  the  integration  step-size.  We  note  that 
Lie-Poisson  systems  appear  naturally  in  many  areas  of  physical  science  and  engineering, 
including  theoretical  mechanics  of  fluids  and  plasmas,  satellite  dynamics,  and  polarization 
dynamics.  In  the  present  paper  we  consider  a  series  of  progressively  complicated  examples 
related  to  rigid  body  systems.  We  also  consider  a  dissipative  example  associated  to  a  Lie- 
Poisson  system.  The  behavior  of  the  mid-point  rule  and  an  associated  reconstruction  rule 
is  numerically  explored. 
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Almost  Poisson  Integration 
of  Rigid  Body  Systems 


1  Introduction 

Natural  dynamical  systems  often  display  a  variety  of  analytic  and  geometric  structures 
in  their  mathematical  descriptions.  Associated  to  such  structures,  there  are  various  con¬ 
served  quantities  or  invariants  of  analytic  as  well  as  geometric  character.  For  instance,  in 
hamiltonian  mechanics  of  particles  in  a  central  force  field,  one  has  conservation  of  energy 
and  total  angular  momentum.  The  phase  space  volume  is  also  conserved.  The  latter  is 
an  example  of  a  geometric  conserved  quantity.  In  the  case  of  dissipative  systems  the  vol¬ 
ume  decays.  In  the  study  of  such  systems  via  computer  simulation,  it  is  very  desirable  to 
use  computational  schemes  that  admit  the  same  set  of  invariants  (or  decay  rates).  More 
precisely,  one  would  like  to  preserve  underlying  geometric  structures  and  symmetries  even 
in  the  discrete-time  dynamics  (computational  scheme),  in  the  interest  of  long-term  pre¬ 
dictions.  If,  as  is  customary,  one  were  to  use  off-the-shelf  schemes  (e.g.  fourth  or  higher 
order  explicit  Runge-Kutta  [21],  backward  Euler,  diagonally  implicit  Runge-Kutta  [2])  to 
integrate  the  dynamics  in  such  problems,  then  the  computed  trajectories  show  systematic 
deviations  (decay  or  growth)  in  the  quantities  that  are  physically  conserved.  Thus,  such 
numerical  simulations  are  an  unreliable  guide  to  the  long-term  dynamic  behavior.  For 
related  comparisons,  see  Channel  and  Scovel  [6]. 

For  some  time,  there  has  been  steady  interest  in  the  design  of  algorithms  that  have  the 
facility  to  closely  mimic  hamiltonian  dynamics.  In  the  work  of  the  Beijing  school  [10,11],  we 
find  a  systematic  exploration  of  symplectic  schemes  (via  generating  functions)  for  classical 
hamiltonian  systems  on  flat  spaces.  In  particular,  the  mid-point  rule  plays  a  prominent  role 
in  that  work,  as  explained  below.  In  the  present  paper,  we  are  concerned  with  systems  that 
evolve  on  cotangent  bundles  of  Lie  groups,  a  basic  example  being  the  free  rigid  body  with 
the  rotation  group  as  configuration  space.  If  the  hamiltonian  is  fully  reducible  by  the  group, 
as  is  the  case  in  the  rigid  body  example,  then  the  dynamics  drops  to  the  flat  space  of  linear 
functionals  on  the  Lie  algebra  of  the  Lie  Group,  and  hence  the  mid-point  rule  is  well-defined 
globally  in  the  reduced  variables.  The  hamiltonian  structure  of  the  reduced  equations  is 
however  noncanonical  and  is  referred  to  as  a  Lie-Poisson  structure.  It  is  a  principal  goal  of 
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this  paper  to  explore  the  applications  of  the  mid-point  rule  and  related  reconstruction  rules 
to  such  Lie-Poisson  equations  via  a  series  of  progressively  complicated  examples.  In  general, 
the  mid-point  rule  is  not  exactly  Poisson  structure  preserving  (and  hence  the  terminology 
of  almost  Poisson  integration).  However,  by  a  small  miracle  involving  the  Jacobi  identity, 
the  mid-point  rule  is  indeed  second  order  accurate  in  the  Poisson  structure  and  to  this  end 
we  give  an  error  formula. 

The  structure  of  the  paper  is  as  follows.  In  Section  2  we  present  the  basic  model  of 
(noncanonical)  Poisson  dynamics  and  summarize  some  of  its  properties.  We  also  specialize 
it  to  the  classical  canonical  setting.  In  Section  3  we  discuss  the  mid-point  rule  and  present 
an  error  formula  for  the  Poisson  bracket.  The  mid-point  rule  is  applied  to  a  variety  of 
examples  in  Section  4;  these  include  rigid  body  dynamics,  heavy  top,  dual  spin  problem, 
and  dual  spin  with  damping.  A  key  result  is  the  derivation  of  a  reconstruction  formula  for 
elements  in  50(3)  (the  rotation  group),  which  conserves  spatial  angular  momenta.  Sections 
5  and  6  discuss  issues  in  the  numerical  implementation  of  the  proposed  algorithms.  Finally, 
numerical  examples  and  simulation  as  well  as  animation  results  for  rigid  body  systems  are 
presented  in  Section  7. 

We  note  in  the  previous  work  of  Marsden  and  Ge-Zhong  [13]  a  Lie-Poisson  Hamilton- 
Jacobi  theory  has  been  developed.  This  theory  leads  to  algorithms  that  preserve  the  Lie- 
Poisson  structure  exactly,  but  do  not  typically  conserve  the  Hamiltonian.  Recently,  Simo 
and  Wong  [24]  have  used  a  Newmark-Based  algorithm  to  study  rigid  body  dynamics.  When 
the  parameters  of  their  algorithms  are  set  so  that  energy  and  momentum  will  be  conserved, 
their  algorithm  reduces  to  the  mid-point  rule  with  an  exponential  map.  We  note,  however, 
that  the  work  of  Simo  and  Wong  does  not  consider  errors  from  the  viewpoint  of  conserving 
the  Poisson  structure,  and  has  not  been  extended  to  applications  with  a  general  Lie-Poisson 
setting.  For  other  earlier  work  on  the  midpoint  rule  we  refer  the  reader  to  Elliot  [8,9]. 

2  The  Model 

In  the  present  paper  we  are  concerned  with  hamiltonian  models  of  the  form, 

z  =  A(z)VH(z).  (1) 

Here  A(z),  the  Poisson  tensor,  is  an  n  x  n  skew- symmetric  matrix  for  each  ^  and  H  is  the 
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hamiltonian. 

identity): 


In  addition,  the  tensor  A(z)  satisfies  a  set  of  differential  equations  (Jacobi 


L. 

i 


dzi 


UZl 


(2) 


With  condition  (2)  the  operation, 

{f,cj}  =  X7fTA(z)Vg 


(3) 


is  a  well-defined  Poisson  bracket  on  the  space  of  smooth  functions  on  IR”.  It  is  remarkable 
that  the  equations  of  many  physical  and  engineering  problems  are  based  on  models  of  the 
form  of  (1)  or  perturbations  thereof.  We  note,  for  example,  the  dynamics  of  dual  spin 
satellites  [17],  accelerator  dynamics  [7],  motion  of  a  heavy  top  [3],  Euler’s  elastica  [20],  and 
dynamics  of  rods,  plates  and  shells  [23].  Also  see  Simo,  Marsden  and  Krishnaprasad  [23] 
for  infinite-dimensional  examples,  and  Marsden  et  al.  [19]  for  a  general  discussion. 

The  dynamics  (1)  can  be  re-expressecl  in  terms  of  the  Poisson  bracket  (3)  as 


=  Ui,  H}.  (4) 

When  A(z)  is  linear  in  z,  the  bracket  structure  is  said  to  be  of  Lie-Poisson  type.  Several 
of  the  earlier  mentioned  examples  are  of  this  variety.  In  this  case  one  sets, 

n 

A,J(d^r^  (5) 

fc-i 

where  T-  =  —Pj,  and  the  Jacobi  identity  (2)  takes  the  form, 

n 

Erfrit  +  r'tri'.  +  r;,r;y  =  o,  i  <  ij.k  <  „.  (6) 

/=1 

It  is  further  well-known  that  in  this  case  the  underlying  vector  space  IR"  can  be  given  the 
structure  of  a  Lie  algebra  [27]  with  the  structure  constants  PP  in  a  suitable  basis.  Now  let 
$  :  lRrt  — »  ]Rre  be  a  diffeomorphism.  $  is  said  to  be  a  Poisson  automorphism  if  it  preserves 
the  Poisson  structure,  i.e.  for  smooth  functions  /,  g , 

{/,</}  o  $  =  {  /  o  $,  g  o  $  }, 
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or  alternatively  its  Frechet  derivative  satisfies, 


(D$(z))A{z)(D$(z)f  =  A($0)). 


(7) 


Here  the  superscript  T  denotes  matrix  transpose.  We  recall  the  well-known  result(see 
Weinstein  [27]), 

Proposition  1  :  The  flow  of  the  system  (1)  satisfies: 

(a)  It  is  a  Poisson  automorphism  V  t. 

(b)  H(&(z))  =  H(z). 

(c)  If  C  :  IR."  — >  1R  is  a  function  such  that 

A  (z)VC(z)  =  0, 

then  C  yields  a  kinematic  conservation  law  of  (1). 

Remark  1:  Functions  C  as  in  Proposition  (1)  are  called  Casimir  functions.  For  canonical 
hamiltonian  systems, 

.  [0  l' 

~  [_i  o  ' 

Here  1  denotes  the  identity  transformation  on  HI”1,  where  m  =  n/2.  The  corresponding 
canonical  Poisson  bracket  is  given  by  the  well-known  formula  [12] 


’  =  E 

i=  1 

where  the  position  coordinates  q,  and  conjagate  momenta  pt  are  given  by 


Of  dg  df  dg 

.  dq,  dpi  dpi  dqi 


z  =  (</!?•••,  <Zrn,  Pi-,--  ■  ,Pm )T  with  n  =  2m. 

For  canonical  hamiltonian  systems,  Casimir  functions  are  constants.  However,  in  the  Lie- 
Poisson  setting,  nontrivial  Casimir  functions  are  common,  as  will  be  demonstrated  by  our 
examples  in  Section  4. 
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3  The  Mid-point  Rule 


A  basic  concern  of  this  paper  is  to  investigate  numerical  algorithms  that  closely  mimic 
Proposition  (1).  In  particular,  we  are  interested  in  the  mid-point  rule,  a  scheme  well-known 
to  be  symplectic  in  the  canonical  case  [11],  Consider  the  implicit  recursion, 

.  h  _ 

This  is  a  discrete  analog  of  (1)  for  time-step  h.  It  is  a  second  order  accurate  integrator, 
and  for  small  enough  h,  defines  a  diffeomorphism,  $hH  via, 

zk+i  =  (zk)  .  (10) 

We  compute  the  Frechet  derivative  D{ ]>hH  (z)  as  follows.  By  definition,  y  =  (c)  is  the 

unique  solution  to  the  implicit  equation, 

F(s,y)=y  -  *  -  h,\(ZZtJL)VH 

=  0.  (11) 


Differentiating  F  (2,  (r))  =  0,  we  get, 

DXF  +  D2  F  o  D  =  0,  (12) 

where  D{F,  i  =  1,2  denote  the  partial  Frechet  derivatives.  For  h  small  enough,  D>  F  has 
an  inverse,  and  (12)  may  be  rearranged  to  give, 


D$hH  (z)  =  -  (D2F)-1  ■  [D1F). 


For  the  special  case  A(z)  =  A  a  constant,  it  is  easy  to  see  that, 


Dif= 


and  DoF  —  1  -  —AH 


+  n  (*) 


(13) 


Letting  Q(z)  =  Hzz 


^  denote 


the  symmetric  Hessian  matrix,  we  get 


/  *  — ^  '  7 

Di>$,  (z)  =  1  -  Q(z)  1  +  -AQ(z)  .  (14) 

Proposition  2:  (Wang  Dao  Liu  [26]).  If  A(z)  =  A  a  constant,  then  the  mid-point  rule 
is  a  Poisson  automorphism. 

Proof:  We  need  to  show  that, 


Di>hH  (z)A(D$hH  (z))T  =  A. 

Based  on  the  above  calculations,  this  reduces  to  showing  that, 

r  _  r  T 

(1  -  \hQ(z))-'  (1  +  \_KQ(z))  A  (1  -  |a Q(=)r'  (1  +  f  A<3(.-)) 

=  A. 

which  is  equivalent  to  showing  that 

h  -]  r  /  -I  T  r  /  If/  1(T 

1  +  -  A  Q(z)  A  1  +  i  A  Q(z)  =  1  -  |  A  Q(z)  A  1  -  ^  A  Q(z)  . 

This  follows  from  the  fact  that  A  =  ~Ar  and  Q(z  )  =  Q(  z)T .  The  proof  is  now  complete. 

Remark  2:  It  follows  from  the  above  proposition  that,  if  A(z)  =  A  =  ^  ^  Q^’i'e'  we 

are  in  the  canonical  case,  then  the  mid-point  rule  preserves  the  classical  Poisson  bracket  (8). 
Thus  we  recover  the  well-known  result  that  the  mid-point  rule  is  a  symplectic  integrator 
(Feng  Kang  [11]).  This  result  also  follows  from  the  observation  that  in  the  canonical  case, 
formula  (14)  becomes  a  Cayley  transform  of  an  infinitesimally  symplectic  matrix  A Q(z) 
and  hence  D$hH{z)  is  symplectic. 

When  A(z)  is  not  a  constant,  in  general,  the  mid-point  rule  is  not  a  Poisson  automor¬ 
phism.  However,  if  A(s)  is  linear  in  z  -  i.e.  we  are  in  the  Lie-Poisson  setting  -  then  the 
following  theorem  shows  that  the  mid-point  rule  is  an  almost  Poisson  integrator  in  the 
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sense  that  it  preserves  the  Lie-Poisson  structure  up  to  second  order.  To  this  end,  an  error 
formula  is  given. 

Theorem  1:  For  Lie-Poisson  systems  the  mid-point  rule  is  almost  a  Poisson  automorphism. 
We  have  the  error  formula, 


D$hH(z)A(z)(D$hH(z))T 

h3 


K(z)A  A( 


A(*k(*)) 

z  +  $hH(z). 


)w)K(zf  +  0(h 4) 


(15) 


where 


A»  =  A( '  +  )Q  + 

Q  =  D':H( 


Here  Q  is  the  Hessian  of  H  evaluated  at  the  mid-point,  and  Q  is  defined  by  requiring  that 
£l(v)y  =  A(y)u,  where  v,  ye  It". 

Proof:  First  note  that,  from  the  definitions  for  partial  Frechet  derivatives,  for  ve  1R." , 


DiF(z,  y)v  A 
D2F{z ,  y)v  A 


lim  F(z  +  tv,  y)  -  F  {z,  y) 

no  / 

lim  F(z,  y  +  tv)  -  F{z,  y) 
t  i  0  t 


and  it  can  be  verified  that, 
DXF=  - 


1  + 


and 


Fb-F 


1 


/l 


+  w 


A(- — &  H(—L 


+  y 


'l 


(10) 

(17) 
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In  the  following  derivation,  y  in  (16)  and  (17)  will  be  used  to  denote  $hH(z).  The  error  in 
the  Poisson  structure  A  due  to  discrete  time  stepping  by  is  given  by, 


£  =  mhu  (z)  A  (z)  (D<b''H(z))T  -  A  («{,(*)). 

Substituting  (13),  (16),  and  (17)  gives, 

£  =  (D2  F)-1  (ei)  (D2  F)t  -1  (19) 

where, 

£!  =  (Dx  F)  A  (z)  ( D !  F)r  -  (D2  F)  A  (y)  (D2  F)T .  (20) 

Since  ( D2  F)~l  =  1  +  0(h),  we  are  mainly  interested  in  the  behavior  of  s i  as  a.  function 
of  h.  For  convenience,  let  us  denote  the  midpoint  {z  +  (*))/ 2  as  x.  In  formula 

(20),  multiplying  out  the  terms,  invoking  (16),  (17)  and  the  linearity  of  A  and  in  their 
respective  arguments,  we  get, 

ti  =  [—  A(<5)  +  hQ(u’) A(,r)  +  h  A(«)fl(tt>)r] 

+  [hA(x)Q  A  +  h  A (x)Q  A(.t)t] 

h2  h 2 

+  -  —  A (x)Q  A(S)QA(x)t  -  jFI(w)A(8)QA(x)t 

h2  h 2 

—  —  A(x)QA(8)Q(w)T  +  —Q(iv)A  (6)Q(w)t  (21) 

Here  8  =  y  —  z  =  (z)  —  The  second  square  bracket  in  t j  vanishes  because 

A  =  —  A7.  Recall  from  (9)  that, 

8  —  liA(x)uj.  (22) 

Substituting  from  (22)  in  (21)  the  first  square  bracket  in  (21)  is  linear  in  h,  and  when 
multiplied  by  v  takes  the  form, 
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h  —  A(A (x)w)  +  £l(w)  A(x)  +  A (x)Q(w)T  v. 


(23) 


We  wish  to  show  that  the  expression  (23)  also  vanishes.  Note  that  this  is  equivalent  to 
showing  that 


[— A(A(ar)iw)  +  fl(w)  A(x)  +  A(a;)f2(u>)r] v  =  0,  V  v  e  IRC  (24) 

Recalling  that  the  definition  for  0  is  given  by  Cl(v)y  =  A(y)v,  the  left  hand  side  of  (24) 
can  be  re-written  as, 


—  Q(v) A  (x)w  +  f l(w)  Q(v)x  +  A(x)  Q,(w)tv.  (25) 

Expanding  0  and  A  via  the  structure  constants  T-),  one  can  show  that  precisely  because 
of  the  Jacobi  identity  (6),  the  expression  (25)  vanishes  for  all  u,  u>,  x  eld".  This  is  the 
small  miracle  alluded  to  in  the  Introduction!  Collecting  together  the  terms  in  the  third 
square  bracket  in  formula  (21)  for  ci  we  get, 

£i 


It  follows  that 

e 


This  completes  the  proof. 

Conserved  Quantities  of  Poisson  System  (1):  It  is  well-known  fact  that  mid-point 
rule  is  a  second  order  algorithm.  Accordingly,  a  conserved  quantity  is  expected  to  be 


-£■  [A(a;)  Q  +  Cl{w)]  A  (A(x)w)  [A(ar)  Q  +  Q(iv)]T 


h3 


K{z )  A 


mL±KM)w 


K(z) 


T 


{D2  F r1  (ei)  (D2  Ff-1 

(1  +  0(h) r1  (ti)  (1  +  0(h))-1 

K(z)a(a(Z  +  ^Z))ro]  K(zf  +  0(h 4) 
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approximated  accurately  up  to  second  order.  In  fact,  letting  F  be  a  first  integral  of  (1),  we 
have 

VF(z)TA(z)VH{z)  =  0,  Vzeir.  (26.1) 

Assuming  F  is  three  times  continuously  differentiable,  by  Taylor’s  theorem  we  can  expand 
F  around  zk  as 


F(zk+1)  =  F(zk)  +  VF(zk)T(zk+1  -zk)  +  -D2F{zk)  •  ( 


k  \  (  ^ k -}- 1  *,k  \  (  „ k -{- 1  _ 

k-\-l 


+  i D3F(zk )  ■  ( zk+1  -  zk)  ■  (zk+1  -  zk)  ■  ( zk+1  -  zk)  +  0(\\zk 
Substituting  the  mid-point  rule  (9)  in  (25.3),  it  can  be  checked  that 


b) 

y  ^  ||^) 

(26.2) 


F(zk+1)-F(zk)  =  ±-D3F(zk)-u-u-u  +  0(h4), 

z4 


(26.3) 


where 


~  k  -|- 1  I  yk-\- 1  I  k 

u  =  AC- . ^  )Vg(— 


)• 


Equation  (26.3)  is  an  error  formula  for  conserved  quantities  of  (1),  which  contains  only 
third  or  higher  order  terms.  It  follows  that  the  mid-point  rule  (9)  preserves  exactly  any 
conserved  quantity  having  only  linear  and  quadratic  terms,  including  Casimir  functions 
and  the  Hamiltonian  for  (1).  As  a  consequence,  we  have  the  following  Proposition  which 
will  be  referred  to  later. 


Proposition  3:  The  discrete  analog  (9)  conserves  all  Casimir  functions  and  the  Hamilto¬ 
nian  H  of  (1)  if  they  contain  only  linear  and  quadratic  terms. 

In  the  next  section,  we  consider  a  set  of  examples  related  to  rigid  body  mechanics  that 
captures  the  essence  of  the  results  of  this  section. 
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4.  Examples  from  Rigid  Body  Mechanics 

The  noncanonical  hamiltonian  model  of  this  paper  -  in  particular  the  Lie-Poisson  case  - 
arises  naturally  in  a  variety  of  problems  in  rigid  body  mechanics.  In  this  section,  we  discuss: 
(a)  the  simple  rigid  body,  (b)  the  heavy  top,  (c)  dual  spin  satellite,  and  (d)  dual-spin  with 
damping. 

4.1  Simple  Rigid  Body  Spinning  Freely  in  Space 

Recall  that  the  equations  of  motion  for  a  simple  rigid  body  in  three  dimensions  take 
the  form: 


A  —  AQ  where  A  £  SO(3)  (27) 

IQ  =  IQ  x  Q.  (28) 


Here,  A  £  SO( 3)  the  group  of  3  x  3  rotation  matrices,  is  the  matrix  of  direction  cosines  for 
a  body  frame  attached  to  the  rigid  body  as  viewed  in  an  inertial  frame;  see  Figure  [1].  The 
vector  Q,  is  the  body  angular  velocity  of  the  rigid  body  and  Q  represents  a  skew- symmetric 
matrix: 


ft  = 


0  -ft3 
ft3  o 
— ft  2  ftj 


-ftl 

0 


The  symbol  x  denotes  the  cross-product  in  H3.  Notice  that  for  an  arbitrary  vector  y, 
f ly  =  0,  x  y.  The  matrix  /  denotes  the  moment  of  inertia  tensor  of  the  rigid  body  in  the 
body  frame.  If  in  =  IQ  denotes  the  body  angular  momemtum,  then  we  can  re-write  (27) 
and  (28)  as: 


A  =  AI~xm  (29) 

777.  =  (30) 


Several  observations  are  in  order: 

[a]  Equation  (30)  is  the  reduction  of  the  full  dynamics  (justified  by  the  S()(  3 )dnvariance 
of  the  kinetic  energy),  [1], 
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Figure  [l]  :  Configuration  of  Inertial  and  Body  Frames 

[b]  Equation  (30)  is  Lie-Poisson  in  the  sense  of  Section  2  with  Poisson  tensor  A (m)  =  m 

and  hamiltonian  H(m)  =  ■  I~lm.  Hereafter  •  denotes  the  dot  product. 

[c]  C(m )  =  A. |jr77. j|2  is  a  Casimir  function  associated  to  (29).  It  is  a  quadratic  conserved 
quantity. 

It  follows  that  if  we  can  integrate  (30),  then  the  kinematic  variable  A  can  be  reconstructed 
by  quadrature  of  (29).  This  process  can  be  mimicked  in  discrete  algorithms  also.  In 
particular,  one  might  integrate  (30)  numerically  via  the  mid-point-  rule  (thereby  conserving 
H  and  C ,  see  Proposition  2)  and  then  devise  a  reconstruction  rule  for  updating  the  attitude 
matrix  in  such  a  way  that  spatial  angular  momentum  tt  =  Am  is  also  conserved. 

4.1.1  Discrete  Update  in  Body  Momentum 

A  two-step  process  is  employed  to  compute  the  discrete  update  in  body  momentum 
and  spatial  attitude.  First,  let  mk  and  rrik+\  be  the  body  momentum  at  timesteps  hk  and 
hk+ 1,  respectively,  and  h  =  hk+i—hk-  The  discrete  update  in  body  momentum  corresponds 
to  the  mid-point  rule  applied  to  (30),  i.e.: 

mk+i-TTik  f  mk  -t  mjk+i  ]  w  r_t 
h  [  2 


rnfc  4-  rrik+i ' 
2 
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(31) 


Equations  (31)  are  solved  by  iteration  for  the  discrete  update  in  body  momentum  mk  to 
mk+i.  It  follows  from  Proposition  3  that  this  update  scheme  conserves  both  H{m)  and 
C(m). 

4.1.2  Discrete  Update  in  Spatial  Attitude  (reconstruction  rule) 

Once  equations  (30)  are  solved,  our  reconstruction  rule  for  the  update 
tude  is  given  by  the  explicit  set  of  equations: 

in  spatial  atti- 

Ak+ 1  —  Ak  1  —bk  J+frjt 

(32.1) 

,  ,  !  \mk  +  mk+1' 

where  bk  —  ^  I 

(32.2) 

The  explicit  recursion  (32)  -  Cayley  Transform  for  50(3)  -  is  arrived  at  by 
servation  of  spatial  angular  momentum  7r  =  Am, 

requiring  con- 

Kk+ 1  =  *k  =  Ak+imk+i  =  Akmk. 

(33.1) 

This  suggests  a  suitable  form  for  the  update  is 

Ak+ 1  —  AkAi  =  Ak  l—bk  1 +bk  , 

(33.2) 

where  bk  is  skew-symmetric  and  needs  to  be  determined  so  that  (33.2)  closely  approximates 
the  matrix  exponential  implied  by  A  =  AQ.  Substituting  (33.2)  into  (33.1)  and  rearranging 
terms  gives 

[??H.  +  mjfc+i]&jk  =  [mk+i  -  rnk]. 

(33.3) 

Notice  that  if  [rnk  +  mk+i]bk  =  [??^r-+i  —  ???./;]  has  a  solution,  then,  by  the  Fredholm  alter¬ 
native  theorem, 

T 

(:>>ik+i  -  mk)  T  kernel  [mk+1  +  rnk] 

(33.4) 

—  kernel  [??2fc+i  +  mk] 

(33.5) 

=  {a[mk+1+mk)  |  aeH}. 

(33.6) 
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This  implies  \mk  -f  mfc.pi]-  [mt+i  —  ??rt]  —  0,  and  thus  conservation  of  Casimir  is  a  necessary 
condition  for  conservation  of  spatial  angular  momentum.  The  general  solution  to  (33.3)  is 
given  by 


h  =  A  t  •  [mk+ 1  -  mk]  x  [mk+i  +  mk]  +  ak[mk+i  +  mk]. 


(34) 


How  do  we  choose  Xk  and  ak  ?  First,  note: 

[mk+ 1  -  rnk]  =  [mk  +  mk+l]bk  (35.1) 

=  [mk  +  mk+i]  x  bk  (35.2) 

=  [mk  +  mfc+i]  x  ([mt+i  -  mk]  x  [mt+i  +  mt])  + 

[mk  +  mk+ i]  x  ajfc[m*+i  +  mfc]  (35.3) 

=  hi  +  mfc+i]  x  \k[mk+1  -  mk]  x  [?«t+i  +  mk]  (35.4) 

=  ||"U-  +  mk+i  ||2  ■  At  •  [mt+i  -  mk\.  (35.5) 


This  implies 


At 


1 

||?nt  +  mt+i 


(36) 


To  get  a t  yt  0,  the  discrete  update  for  body  momentum  is  substituted  into  (34).  Rearrang¬ 
ing  terms  and  selecting 


Oik  = 


[mk  +  777-t+i]  •  I  1  [mk  +  mt+i] 
||  rnk  +mk+i\\2 


(37.1) 


gives 


bk 


mk  +  mk+1 

9 


(37.2) 


Remark  3:  Table  [1]  summarizes  the  discrete  equations  for  Rigid  Body.  Notice  that  the 
matrix  product 


1-6* 


1+bk 


=  exp 


hi 


-l 


mk  +  mt+i 

9 


Q(h3 


(38.1) 
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In  fact,  if  I  xm  in  (29)  is  a  constant  vector  across  the  timestep,  then  the  incremental  rule 
for  the  attitude  matrix  becomes 

Ak+ 1  =  Ak  exp {hl~lm).  (38.2) 

Thus,  the  update  rule  (33.2)  can  be  intuitively  thought  of  as  the  composition  of  two  steps. 
First,  the  angular  momentum  in  (29)  is  approximated  by  an  intermediate  value.  Second,  the 
exponential  in  (38.2)  is  approximated  in  terms  of  b k  by  the  formula  (38.1).  We  summarize 
our  energy-momentum  preserving  algorithm  for  a  simple  rigid  body  spinning  freely  in  space 
as  follows. 

Update  for  Body  Momentum: 

mk+ i  -  nik  _  \mk  +  mk+1  ]  T_x  \ mk  +  ?7u-+i 

h  -[  2  J  X  1  [  2 

Update  for  Spatial  Attitude: 

_ i — i r 

Ak+i  —  Ak  1  —hk  l+bk 

,  ,  !l  ]  r-i  [  'ilk  +  mk+1 ' 

Table  [1]  :  Discrete  Equations  for  Rigid  Body 
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4.2  Heavy  Top 


Consider  the  motion  of  a  rigid  body  fixed  at  a  stationary  point  subject  to  the  action 
of  a  uniform  gravitational  field,  as  shown  in  Figure  2. 


Figure  [2]  :  Diagram  for  Heavy  Top 


If  the  body  frame  is  fixed  at  a  point  0,  then  we  define  the  unit  vector  from  0  to  the  center 
of  mass  in  the  body  frame  as  y.  The  spatial  inertial  frame  is  also  at  0  with  the  vector 
k  being  one  of  its  axis.  Notice  that  the  body  is  acted  on  by  a  gravitational  force  —Mg k, 
where  M  is  the  mass  of  the  rigid  body,  and  g  is  the  constant  of  acceleration  due  to  gravity. 

The  kinematics  of  the  body  frame  relative  to  the  interial  frame  are  governed  by 

A  =  An  where  A  £  50(3).  (39) 

Moreover,  if  v  —  AT  ■  k,  then  the  equations  of  motion  for  the  heavy  top  can  be  expressed 
as 


m  =  m  x  I  1m  +  Mglv  x  y 

(40) 

v  —  v  x  I~lm 

(41) 

where  m  =  in  is  body  momentum,  l  the  distance  between  0  and  the  body  center  of  mass, 
and  I  the  moment  of  inertia  of  the  heavy  top  relative  to  the  body  frame.  In  addition,  we 
note: 
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a]  Equations  (40)  and  (41)  are  Lie-Poisson  in  the  sense  of  Section  2  with  Poisson  tensor 


A  (m)  — 


v  0 


and  total  energy 


H  (m,v)  —  —m  ■  I  1m  -f  Mglv  ■  x- 


See  Marsden,  Ratiu  and  Weinstein  [IS]  for  details. 

[b]  The  heavy  top  problem  has  two  Casimirs: 

Cl  (m,  v)  =  i||u||2  (43) 

Co  (m,  v)  —  m  ■  v.  (44) 

Equation  (43)  reflects  the  fact  that  .4  £  50(3),  and  equation  (44)  says  that  spatial 
angular  momentum  along  the  vertical  axis  k  passing  through  O  is  conserved. 


4.2.1  Discrete  Update  for  Heavy  Top 

The  corresponding  discrete  equations  of  motion  using  mid-point  rule  are 

rrik+i  -  mk  f  m*  +  mjt+i  ]  w  T_x  mk  +  mk+i  ,  *,r r.,vk  +  Vk+ 1 

q  ■  A'lgi  x  x 


Vk+l  ~  vk  I  vk  +  vk+x 


in  k  +  riik+i 


(45.1) 


As  noted  in  Proposition  3,  the  mid-point  rule  preserves  exactly  the  Hamiltonian  and 
Casimirs  which  contain  only  linear  and  quadratic  terms.  The  above  scheme  (45)  conserves 
both  Casimir  functions  C i,  C-2  in  (43),  (44)  and  the  Hamiltonian  in  (42). 

4.2.2  Discrete  Update  in  Spatial  Attitude 

After  equations  (45)  are  solved,  the  discrete  update  and  spatial  attitude  of  the  heavy 
top  is  given  by  the  equation 


-4r+i  =  Ak  1— h  1+fr/t 


(46.1) 


where  h=  f-lj-1  [mt  +  mt+1 

0  9 
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The  justification  of  the  reconstruction  rule  (44)  is  analogous  to  the  one  used  for  the  simple 
rigid  body. 

4.3  Dual  Spin  Satellite 

Dual  spin  satellites  consist  of  a  simple  rotating  platform  carrying  one  or  more  sym¬ 
metric  rotors  spinning  at  constant  rates  about  an  axis  fixed  relative  to  the  platform.  The 
purpose  of  the  rotors  is  to  exert  internal  reaction  torques  on  the  platform.  Assuming  there 
are  no  external  forces  acting  on  the  system,  the  net  effect  is  a  transfer  of  interned  momen¬ 
tum  from  the  platform  to  the  rotors.  In  the  presence  of  a  suitable  damping  mechanism 
and  sufficiently  high  rotor  velocities,  the  spacecraft  angular  velocity  will  align  itself  with 
the  rotor  axis. 


Figure  [3]  :  Schematic  of  Dual  Spin  Satellite 


For  the  purposes  of  this  study  we  are  interested  only  in  dual  spin  satellites  where  the  axis 
of  the  rotors  passes  through  the  center  of  mass  of  the  platform  with  fixed  rotors.  The 
attitude  kinematics  of  the  satellite  relative  to  an  inertial  frame  is  given  by 

A  =  AQ.  where  A  €  50(3).  (47) 

The  equation  of  motion  for  body  angular  momentum  is  now, 

m  =  (m  +  0  •  Vmi7(m) 

=  [m  +  l]  x  I_1m 
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(48.1) 

(45.2) 


where  m  is  a  vector  of  body  angular  momentum,  and  I  is  the  locked  moment  of  inertia 
tensor  of  the  body  plus  rotor  system.  The  satellite  has  an  internal  rotor  spinning  at  constant 
relative  angular  momentum  /.  We  also  note: 

[a]  Equations  (48)  are  Lie-Poisson  in  IR3  with  Poisson  tensor  A (m)  =  in  +  /,  and  hamil- 

tonian  H  (m)  =  •  I~l  m.  See  Krishnaprasad  [17]  for  details. 

[b]  The  dual  spin  satellite  has  one  Casimir  C  (m)  =  ||m  +  /||2. 

[c]  The  spatial  angular  momentum  is  ir  —  A  {in  +  /),  and  is  conserved  since  the  torques 
generated  by  the  rotors  are  internal  torques. 


4.3.1  Discrete  Update  in  Body  Momentum 

The  discrete  update  in  body  momentum  is  given  by  the  midpoint  rule: 


mk+i  -  mk 
h 


mk  +  mk+ 1  +  ^ 


x  r 1 


mk  +  mk+i 

9 


(49) 


4.3.2  Discrete  Update  in  Spatial  Attitude 

The  discrete  update  in  spatial  attitude  is 


Ak+ 1  —  Ak 


1  -h 


i  -1  r 


1  +bk 


where 


bk  = 


'  h 

T~l 

mk  +  mk+1 

2 

1 

2 

(50.1) 

(50.2) 


Remark  5:  Following  the  steps  in  Section  4.1  it  is  easy  to  verify  that  the  attitude  update  is 
momentum  conserving.  Indeed,  techniques  for  reorienting  the  attitude  of  a  satellite  depend 
on  a  mechanism  for  transfering  the  balance  of  internal  rotor  and  platform  momenta,  without 
affecting  spatial  angular  momentum.  Independently,  by  Proposition  3,  the  energy  H  and 
Casimir  C  are  conserved  by  the  momentum  update  (49). 

4.4  Dual  Spin  Satellite  with  Damped  Rotors 

Dual  spin  satellites  consist  of  a  rigid  body  platform  and  several  internal  rotors.  When 
the  platform  is  fully  operational,  the  rotors  are  set  to  spin  at  a  constant  angular  velocity 
relative  to  the  platform.  Damping  rotors  act  as  dissipators  of  energy.  Indeed,  even  in  the 
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Figure  [4]  :  Schematic  of  Dual  Spin  Satellite  with  Damping 


presence  of  mild  disturbances,  the  transfer  of  momentum  from  the  platform  to  the  damping 
rotors  will  result  in  a  mechanism  for  attitude  aquisition  of  the  platform. 

The  equations  of  motion  for  a  dual  spin  satellite  with  damping  are  composed  of  two  parts. 
As  with  the  two  previous  applications,  the  kinematics  are  given  by 

A  =  Afi!  where  A  €  SO(3).  (51) 

Let  the  diagonal  matrix  a  =  diag  (01,0:2,0:3),  where  aq  >  0,  and  let  Id  =  diag  l|) 

be  a  diagonal  matrix  of  moments  of  inertia  of  the  damping  rotors  with  respect  to  their  spin 
axes.  Note  that  /^  >  0  for  i  =  1, 2,  3.  Moreover,  if  m  and  l  are  as  previously  defined,  and 
d  a  momentum  vector  associated  to  damping  rotors,  then  the  evolution  of  body  angular 
momentum  is  governed  by  the  equations 

rh  —  (m  +  l  +  d)  x  I~lm  —  7 m  +  8cl  (52.1) 

d  =  7  m  —  8d  (52.2) 

where  7  =  al~l  and  8  =  alf1,  respectively.  Although  the  system  is  not  a  hamiltonian 
system  in  the  form  of  (1),  we  still  have  the  following  properties  (see  Krishnaprasad  [17]  for 
details): 

[a]  This  system  has  one  conserved  quantity  C  (m,d)  =  || m  +  l  +  d||2. 
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(53.1) 


[b]  The  Lyapunov  function  for  dual  spin  with  damping  is: 

V  ( m,d )  =  —m  ■  I~1?n  +  —cl  ■  Ici~1d 


It  is  decrescent  with 


V 


(Id~1d  -  rlm)r  a(Ici~1d  -  r'm). 


1~\T  tT  — 1 


(53.2) 


along  trajectories  of  (52).  Again,  see  Ivrishnaprasad  [17]  for  details. 

[c]  The  spatial  angular  momentum  tt  =  A  [ m  +  l  +  d]  is  also  a  conserved  quantity. 

4.4.1  Discrete  Update  in  Body  Momentum 

The  update  in  body  momentum  is  given  by  the  mid-point  rule 


rafc+i  ~  mk 
li 


f?U'  +  l  ,  dk  +  dk+ 1 

—  r  t  i  "7 


x  I 


-i 


7 


dk+ 1  -  dk 
h 


=  7 


'mk  +  mjt+i' 

+  6 

dk  +  dk- )-i 

2 

2 

mk  +  mk+ 1 

-  8 

dk  +  dk- i-i 

2 

2 

mk  +  mk+i 


(54.1) 

(54.2) 


4.4.2  Discrete  Update  in  Spatial  Attitude 

The  discrete  update  for  spatial  attitude  is: 


where  hi.  = 


Ak+i  =  A  * 
h 


1— 


1 +bk 

mk  +  mk+ 1 


(55.1) 

(55.2) 


4.4.3  Conservation  Law  for  Dual  Spin  with  Damping 

The  system  (52)  admits  C  and  ir  as  conserved  quantities.  The  Lyapunov  function 
V  decreases  at  a  rate  quadratic  in  the  damping  vector  cl.  The  discrete  scheme  (54)-(55) 
mimics  these  features,  as  one  might  expect.  We  leave  it  to  the  reader  to  verify  that  the 
discrete  version  of  (53.2)  is  satisfied. 
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5.1  Numerical  Integration  Schemes 

The  solution  procedure  for  each  of  the  applications  described  in  Sections  (4.1)  through 
to  (4.4)  is:  (a)  solve  the  implicit  equations  for  the  body  momentum  at  timestep  tk+i ,  (b) 
solve  the  explicit  equations  for  the  update  in  spatial  attitude. 

5.1.1  Solution  of  Implicit  Equations  for  Update  in  Body  Momentum 

Equations  (31),  (45),  (49)  and  (54)  are  a  systems  of  nonlinear  ordinary  differential 
equations  (in  vector  form)  that  must  be  solved  at  each  timestep  for  . 

First,  each  equation  is  written  in  component  form  and  rearranged  so  that  solving  the 
problem  is  equivalent  to  finding  the  root  of  an  equation.  This  gives  6  equations  for  dual 
spin  with  damping  and  heavy  top  applications,  and  3  equations  for  rigid  body  and  dual 
spin.  Equations  (58)  in  Table  [2]  show,  for  example,  the  3  component  equations  of  body 
momentum  for  the  dual  spin  problem.  The  adopted  notation  for  ??^,  l  or  I  is  timestep  i  of 
component  j  in  the  [f,j]  subscript  brackets. 

A  damped  Newton-Raphson  procedure  is  used  to  solve  the  component  equations  at 
each  timestep.  If  is  the  pth  iterate  of  body  momentum  at  timestep  k+ 1,  then  the  equa¬ 
tions  to  be  solved  at  each  iterate  are  obtained  by  expanding  m£+1),  giimi,--  mpk+1 ) 

and  g^irnki  mjj!+1)  in  a  Taylor  series  about  m£+1,  truncating  all  second  order  terms  and 
higher  terms,  and  solving  the  set  of  equations: 


dg  i 

9.91 

dcji 

9.92 

SmF*  + 1,2] 

9.92 

992 

dm[k+ i.i] 

dgs 

dmlk  +  l,2] 

9.9a 

amft+113] 

9.9a 

9mp  +  i.2] 

+  1,3] 

Amfl-,2] 

Amf*,3] 


'  gi(mk,mpk+1)' 
g2(mk,mpk+1) 
_g:Arnk,mpk+1)  _ 


(56) 


for  the  incremental  update: 


<+1.0  =  +  Amffc,ij  (57) 

where  i  =  1,2,3.  Equations  (56)  and  (57)  may  be  stated  more  concisely  J.h  =  — g. 
Iterations  continue  at  each  timestep  until:  (a.)  a  preset  maximum  number  of  iterations 
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Update  for  Body  Momentum: 


9i  (”**,”**+ 1)  =  mp.-+i,i]  -  j 

'h 

(”*[*,  2]  +  ”*[*+1,2]  +  2 

_4_ 

'h 

(”*[*,3]  +  ”*[*+1,3]  +  2 

+ 

_  4_ 

-  0 

g2  (mk,mk+i)  =  m[fc+1|2]  -  m[*,2] 

'h 

(”*[*,  3]  +  ”*[*+1,3]  +  2 

— 

_4_ 

'll' 

(”*[*, i]  "b  J”[*+i,i]  +2 

+ 

_  4_ 

(”*[*,3]  +  ”*[*+1,3]. 


(58.1 ; 


(”*[*, 1]  +  ”*[*+1,1]. 


(”*[*,  3]  +  ”*[*+ 1,3]. 


=  0 

g3  (mk,mk+1)  =  m[k+  tf3j  -  m[k^ 

h 


(58.2) 


(”*[*, i]  +  ”*[*+!, i]  +  2  •  h)  ■ 


(m[k,  2}  +  rn[k+ 1,2 


+ 

=  0 


/  ,  ,  o  ,  \  (m[M]  +  m[*+U]] 

(”*[*.2]  +  ”*[*+1.2]  +  2  ’  h) - 7 - 

J- 1 


(58.3) 


Table  [2]  :  Dual  Spin  Equations  in  Rearranged  Component  Form 


is  reached,  or  (b)  all  changes  in  body  momentum  components  from  (57)  are  less  than  a 
preset  error  value  times  the  magnitude  of  momentum  component  at  the  beginning  of  the 
iterations.  Moreover,  divergence  of  the  iterates  is  avoided  by  ensuring  \\g  (mk,  ”*fc+i)  ||  is 
less  than  \\g(rnk,  mPk+1)  ||2-  When  this  test  fails,  li  is  divided  by  powers  of  2  until  the  norm 
of  the  new  [p  +  lf,!  iterate  is  less  than  ||</(???.fc,  niPk+1 )  ||9- 
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5.1.2  Solution  of  Explicit  Equations  for  Update  in  Attitude 

Numerical  solutions  to  the  equations  in  Table  [1]  are  found  by  first  solving  (31)  im¬ 
plicitly  for  mjt+i.  From  a  theoretical  standpoint,  m^+i  can  be  explicitly  substituted  into 
(32.2),  and  then  (32.1)  for  Ak+ i-  Because  it  is  undesirable  to  compute  the  matrix  inverses, 
solutions  to  (32.1)  are  obtained  by  first  solving 


1-h 


1  +h 


(59.1) 


for  the  (3  x  3)  matrix  [x*],  followed  by  the  explicit  update 

Ak- h  =  Ak  ■  [.Tfc]. 


(59.2) 


6.  Implementation 

The  implementation  of  algorithms  described  in  Sections  4. 1-4.4  has  been  underway  at 
the  Systems  Research  Center  since  19SS  [5,22].  Initially,  these  simulations  were  written  on 
a  Silicon  Graphics  IRIS  3130  workstation  with  integrated  graphical  and  numerical  code. 
While  this  approach  is  satisfactory  for  applications  that  are  computationally  non-intensive, 
Sela  [22]  reports  that  increased  computational  power  would  be  needed  for  computationally 
intensive  applications  to  achieve  simulation  in  real  time.  Since  the  IRIS  Workstation  has 
customized  hardware  for  the  Graphics,  a  natural  solution  was  to  develop  a.  graphically 
based  user  interface  on  the  IRIS  for  the  animation,  and  setting  up  the  problem  description, 
and  run  the  numerical  simulation  as  a  separate  server  process  on  a  fast  SUN  workstation. 
InterProcess  Communication  Facilities  [16]  are  used  to  connect  the  processes  via  Ethernet. 

Design  of  the  user  interface  was  based  on  a  toolkit  developed  by  NAS  A- AMES  [25] 
specifically  for  the  Silicon  Graphics  IRIS  workstation  family.  The  library  consists  of  dials, 
sliders,  and  other  mouse-sensitive  actuators,  and  panels,  which  are  groups  of  actuators  that 
appear  as  separate  windows  on  the  IRIS  workstation. 

When  the  program  is  first  activated,  the  user  interface  process  located  on  the  IRIS  es¬ 
tablishes  an  IPC  connection  with  the  numerical  simulation  server,  sets  up  a  variety  of  panel 
windows  for  the  user  interface,  and  waits  for  the  user  to  initialize  simulation  parameters 
before  sending  the  simulation  request  via  IPC  to  the  numerical  simulator. 
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Figure  [5]  :  User  Interface  for  Rigid  Body  Simulation 


Slider  Name 

Description 

Toggles 

Buttons  to  remove  sliders  from  display. 

Rotor 

Set  3  components  of  angular  momentum  of  internal  rotor. 

Omegas 

Set  3  components  of  initial  body  angular  momentum. 

Timestep 

Set  timestep  length  for  Midpoint  Integration.  Minimum  and  maximum  startup 

default  values  are  0.0  and  0.2  seconds,  respectively. 

Inertia 

Dual  Slider  to  set  total  interia  of  the  body,  and  component  of  rotor  inertia  along 

each  axis.  In  all  cases  the  rotor  inertia  must  be  less  than  the  inertia  of  satellite 

with  fixed  rotors. 

Damping 

Set  damping  ratio  for  rotors.  Notice  that  although  equation  (52)  theoretically 
permits  different  damping  ratios  along  each  axis,  here  we  assume  all  components 

of  damping  are  identical. 

Output 

Displays  output  for  time  variation  of  Casimir,  Total  Energy,  and  each  component 

of  body  angular  momentum. 

Table  [3]  :  Description  of  User  Interface  Sliders 
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Figure  [5]  shows  a  typical  screen  layout  for  the  user-interface.  Individual  panels  con¬ 
taining  single  and  groups  of  sliders  are  provided  for:  (a)  setting  the  integration  timestep, 
(b)  level  of  damping  in  the  rotors,  (c)  components  of  the  principal  moments  of  inertia  of 
the  combined  platform  and  rotors,  (cl)  principal  moments  of  inertia  of  the  rotors,  and  (e) 
components  angular  momentum  for  the  rotors.  In  an  effort  to  provide  insight  into  the 
dynamical  behavior  of  dual  spin  satellites,  we  assume  for  the  purposes  of  graphical  display 
that  the  satellite  platform  is  a  rectangular  block  whose  dimensions  change  with  adjustments 
to  the  prinicipal  inertia  sliders.  Thus,  a  riser  becomes  immediately  aware  that  the  moment 
of  inertia  sliders  cannot  be  set  in  an  arbitrary  manner.  This  is  because  a  solid  rigid  body 
must  satisfy  the  triangle  inequalities 

h  <  h  T  I 2 i  I*<h  +  Iu  and  h  <  h  +  h •  (60) 

The  interested  reader  is  referred  to  Chapter  0  of  Arnold  [3]  for  details.  Also  displayed  are 
the  body  axes  for  the  rectangular  block,  and  the  resultant  momentum  vector  for  the  internal 
rotors.  Users  quickly  determine  that  adjustments  to  the  matrix  of  damping  coefficients,  a, 
have  no  affect  on  the  Casimir  C  (m,  d)  —  jjm  +  l  +  d ||2,  but  do  change  the  rate  of  decay  for 
equation  (53.1).  Conversly,  changes  to  the  angular  momentum  components  for  the  rotors 
result  in  an  almost  immediate  adjustment  of  the  Casimir,  with  no  apparent  jump  in  the 
Lyapunov  function  decay  rate. 

7.  Numerical  Experiments 

Results  of  numerical  experiments  are  presented  for  two  applications:  (a)  motion  of  a 
rigid  body  spinning  freely  in  space,  and  (b)  dynamics  of  a  dual  spin  satellite  platform  con¬ 
taining  high  spinning  internally  damped  rotors.  In  each  case  we  are  interested  in  verifying 
that  simulations  of  the  discrete  dynamics  match  the  theoretical  predictions  of  Section  4. 

7.1  Example  1  :  Rigid  Body  Motion 

Consider  the  motion  of  a  rigid  body  having  principal  moments  of  inertia  I  = 
diag  (1,  2,  3).  Initially  (at  time  t  =  0  )  the  rigid  body  is  spinning  freely  about  its  inter¬ 
mediate  (unstable)  axis  with  angular  velocity  to  =  [1,10,1].  400  timesteps  at  At  =  0.05 
seconds  are  computed.  At  each  timestep,  iterations  to  solve  the  equations  of  motion  for 
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body  momentum  rrik+i  continue  until  the  magnitude  of  the  incremental  updates  in  angular 
momentum  components,  given  by  equation  (57),  are  all  less  than  1/100  of  the  component 
magnitudes  for  load  vector  g  at  the  beginning  of  the  iterations. 

Energy  and  Casimir  Conservation.  Figure  6  shows  the  time  variation  in  Energy 
H(mk+i)  =  \n~ik 4.1  •  and  Casimir  C (»??.*+ 1)  =  \  ||???.£._|_i  ||2  for  the  rigid  body. 

Both  quantities  are  conserved  to  machine  precision. 

Spatial  Attitute.  Figures  7  and  8  show  the  time  variation  in  attitude  matrix  component 
4(1,1)  for  400  steps  of  simulation  at  At  =  0.05  seconds,  and  4000  steps  of  simulation  at 
At  =  0.005  seconds,  respectively.  We  recall  that  the  role  of  the  attitude  matrix  A  is  to 
describe  the  rotational  orientation  of  the  body  frame  relative  to  an  inertial  frame.  As  such, 
the  time  variation  in  all  components  of  A  is  bounded  by  the  interval  [—1,1].  Figures  7 
and  8  match  this  observation.  Our  initial  numerical  experiments  were  with  timestep  At  = 
0.05.  Although  the  Casimir  and  Energy  were  conserved  to  machine  precision,  we  suspected 
the  validity  of  the  sharp  points  in  the  plot  of  .4(1,1)  in  Figure  7.  To  ensure  that  this 
was  not  an  error  in  coding  (it  happens!)  the  simulations  were  repeated  for  the  same  time 
inteval  using  At  =  0.005  seconds.  The  results  are  plotted  in  Figure  S.  It  is  easy  to  see 
that  time  variations  in  peak  values  of  attitude  component  4(1,1)  of  Figure  8  are  much 
smoother  than  Figure  7.  However,  the  dramatic  contraction  in  the  period  of  the  attitude 
components  by  merely  decreasing  the  step  length  was  not  anticipated  apriori.  This  can  be 
seen  by  counting  the  number  of  cycles  of  4(1, 1)  over  t  E  [0,  20]  seconds  for  Figures  7  and 
8.  The  former  has  approximately  28.5  cycles,  and  Figure  S  approximately  29.3  cycles.  We 
note  that  this  observation  is  consistent  with  period  extensions  in  the  Newmark  method;  for 
a  discussion,  see  Chapter  9  of  Hughes  [15].  Unfortunately,  this  observation  also  exposes  the 
main  weakness  of  the  numerical  approximation.  While  Casimirs  and  Energy  are  conserved 
to  machine  precision  (indicating  that  the  discrete  dynamics  follow  the  true  trajectories  in 
reduced  phase  space),  we  know  from  Section  2  that  the  mid-point  rule  is  only  second  order 
accurate  in  computing  the  Poisson  bracket.  This  translates  to  a  systematic  deviation  of  the 
discrete  attitude  from  time-varying  attitude  of  the  continuous  system.  Work  is  currently 
underway  to  try  and  improve  the  attitude  prediction  by  combining  the  mid-point  rule  with 

Richardson’s  Extrapolation  techniques  [4j. 
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7.2  Example  2  :  Dual  Spin  with  Damping 

Our  second  example  examines  the  discrete  time  response  of  a  dual  spin  statellite  plat¬ 
form  containing  3  internal  rotors  spinning  at  constant  angular  velocity  relative  to  the 
satellite  platform.  Principal  moments  of  inertia  for  the  rotors  are  Id  =  dicig  (0.1,  0.1, 0.1). 
After  the  rotors  are  locked,  the  combined  rotors  plus  satellite  platform  is  assumed  to  have 
principal  moments  of  inertia  I  =  cliag  (1, 2, 3).  The  internal  rotor  is  spinning  with  angular 
momentum  components  L  —  [0,0,10].  At  time  t  =  0,  the  rigid  body  is  spinning  freely 
about  its  intermediate  (unstable)  axis  with  angular  velocity  w  —  [1, 10, 1].  2000  timesteps 
at  At  =  0.05  seconds  are  computed. 

Conserved  Quantities.  Our  discrete  approximation  to  the  Lyapunov  function  given  in 
equation  (51)  asymptotically  approaches  a  non-zero  value  as  theoretically  predicted.  At 
the  same  time,  the  Casimir  C  (mk+i ,  <4+1 )  =  ||nu-+i  +  l  +  c4+i||2  is  conserved  to  machine 
precision. 

Spatial  Attitute.  Our  numerical  experiments  conserve  spatial  angular  momentum 
TTk+i  =  A[rrik+ 1  +  /  +  dk+ 1]  to  machine  precision.  As  mentioned  in  Section  4,4,  momen¬ 
tum  transfer  relies  on  the  conservation  of  this  quantity,  and  an  alignment  of  the  platform 
attidude  along  a  single  axis  follows  for  high  enough  speeds  of  the  driven  rotors.  Time 
variations  in  components  of  the  attitide  matrix  .4(2,2)  and  .4(2,  3)  are  shown  in  Figures  10 
and  11,  respectively.  Although  the  platform  is  initially  tumbling  in  an  erratic  manner,  the 
effect  of  damping  results  in  an  alignment  of  the  rotation  about  a  single  axis.  This  is  in¬ 
ferred  from  component  A( 2, 3)  asymptotically  approaching  a  single  value,  while  component 
A( 2, 2)  appears  to  approach  a  steady  oscillatory  motion. 
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Energy  and  Casimir 


Energy  and  Casimir 


Time  (seconds) 


Figure  [8]  :  A(l,l)  vs  Time  with  At  —  0.005  Seconds 


Time  (seconds) 


Figure  [9]  :  Energy  and  Conserved  Quantities  vs  Time  for  Dual  Spin  with  Damping 
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8  Conclusions 


Exact  solutions  to  the  flow  of  a  Lie-Poisson  system  correspond  to  a  continuous  succes¬ 
sion  of  Poisson  automorphisms  [13],  and  conserve  all  Casimirs  and  the  Hamiltonian  (when 
applicable).  This  work  has  been  motivated  by  a  need  to  find  discrete  approximations  to 
the  flow  that  have  the  same  properties.  In  particular,  we  have  focussed  on  the  numerical 
integration  of  Lie-Poisson  systems  using  the  mid-point  rule.  We  have  proved  that  the  mid¬ 
point  rule  preserves  the  Poisson  structure  up  to  second  order  accuracy.  An  explicit  error 
formula  has  been  derived,  and  may  be  used  in  measuring  the  deviations  of  the  Poisson 
structure  through  the  approximated  flow.  Moreover,  we  have  shown  that  any  naturally 
conserved  quantity  can  be  approximated  up  to  second  order  by  using  the  mid-point  rule. 
For  conserved  quantities  containing  only  linear  and  quadratic  terms  such  as  those  in  the 
examples  of  this  paper,  the  associated  conservation  laws  are  exactly  satisfied. 

Four  rigid  body  systems  have  been  studied.  In  each  case,  the  dynamics  in  full  phase 
space  are  obtained  via  an  explicit  reconstruction  rule.  The  result  is  conservation  of  spatial 
angular  momentum  for  the  rigid  body  and  dual  spin  satellites,  and  conservation  of  momen¬ 
tum  along  the  vertical  axis  for  the  heavy  top  example.  Indeed,  the  attitude  acquisition  of 
satellites  containing  damping  rotors  depends  on  an  internal  transfer  of  momentum,  while 
also  conserving  spatial  angular  momentum.  The  behavior  of  dual  spin  satelites  and  rigid 
body  applications  has  been  animated  on  an  IRIS  workstation. 

Future  applications  of  study  will  include  the  approximations  of  the  motion  of  a  satellite 
in  a  central  gravitational  field  and  optimal  control.  Work  is  currently  underway  to  try 
and  improve  the  attitude  prediction  by  combining  the  mid-point  rule  with  Richardson’s 
Extrapolation  techniques  [4],  In  the  long  term,  we  share  the  view  of  Simo  and  Wong  [24] 
that  approximation  precedures  for  the  dynamics  on  SO(3)  is  directly  applicable  to  transient 
dynamic  calculations  of  geometrically  exact  rods  and  shells. 
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